#include "Geometry.h"

using namespace blitz;
void MoM_Geometry::ReadGeometry(std::string GeoFileName, int buffer)// Open Geometry File
{
	std::cout<<"Opening Geometry \""<<GeoFileName<<"\""<<endl;
	fin.open(GeoFileName);									
	if (!fin) std::cerr<<"File not exist!"<<std::endl;
	Buffer_Layers = buffer;
	ImportNode();
	ImportPatch();
	ImportEdge();

	Edge_Max = 0;													// Saving Edge Length
	for (int i=0; i<Edge_Num; i++)
	{
		int node1 = Edge(i)[0];
		int node2 = Edge(i)[1];
		Edge_Length(i) = sqrt(dot((Node(node1)-Node(node2)),(Node(node1)-Node(node2))));
		if (Edge_Length(i)>Edge_Max) Edge_Max = Edge_Length(i);
	}

	for (int i=0; i<Patch_Num; i++)									//Saving Patch Center, Area and Normal Vector 
	{
		int node1 = Patch(i)[0];
		int node2 = Patch(i)[1];
		int node3 = Patch(i)[2];
		Position3D R1,R2,R3;
		R1 = Node(node1);
		R2 = Node(node2);
		R3 = Node(node3);
		Patch_Ctr(i) = (R1+R2+R3)/3;
		Position3D R12,R23,R31;
		R12 = R2-R1;
		R23 = R3-R2;
		R31 = R1-R3;
		Position3D tmp = cross(R12,R23);
		Patch_Norm(i) = tmp/sqrt(dot(tmp,tmp));
		real_data e1 = sqrt(dot(R12,R12));
		real_data e2 = sqrt(dot(R23,R23));
		real_data e3 = sqrt(dot(R31,R31));
		real_data s = 0.5*(e1+e2+e3);
		Patch_Area(i) = sqrt(abs(s*(s-e1)*(s-e2)*(s-e3)));
	}
	GdPtWt.resize(2,4);
	GdPtCoef1.resize(2,4);
	GdPtCoef2.resize(2,4);
	GdPtCoef3.resize(2,4);
	R.resize(2,Patch_Num,4);
	Initial();
	fin.close();
}

void MoM_Geometry::ImportNode()										// Import Node Information
{
	fin>>Node_Num;															
	Node.resize(Node_Num);
	int ii;
	for (int i=0; i<Node_Num; i++)
	{
		fin>>ii>>Node(i)[0]>>Node(i)[1]>>Node(i)[2];
	}
}

void MoM_Geometry::ImportPatch()									// Import Patch Information
{
	fin>>Patch_Num;														
	Patch.resize(Patch_Num);
	Patch_Ctr.resize(Patch_Num);
	Patch_Norm.resize(Patch_Num);
	Patch_Area.resize(Patch_Num);
	int ii;
	for (int i=0; i<Patch_Num; i++)
	{
		fin>>ii>>Patch(i)[0]>>Patch(i)[1]>>Patch(i)[2];
	};
	Patch = Patch - 1;
}

void MoM_Geometry::ImportEdge()										// Import Edge Information
{
	fin>>Edge_Num;
	I.resize(Edge_Num);
	I = 0;
	Edge.resize(Edge_Num);
	Edge_Length.resize(Edge_Num);
	int ii;
	for (int i=0; i<Edge_Num; i++)
	{
		fin>>ii>>Edge(i)[0]>>Edge(i)[1]>>Edge(i)[2]>>Edge(i)[3];
	};
	Edge = Edge - 1;
}

void MoM_Geometry::Initial()
{
	Near2 = 1.02 * Edge_Max * Edge_Max;
	SetCoef();
	for (int Rule=0; Rule<2; Rule++)
		for (int Pat=0; Pat<Patch_Num; Pat++)
		{
			int Node1 = Patch(Pat)[0];
			int Node2 = Patch(Pat)[1];
			int Node3 = Patch(Pat)[2];
			for (int Grid=0; Grid<Grid_Num(Rule); Grid++)
				R(Rule,Pat,Grid) = GdPtCoef1(Rule, Grid) * Node(Node1) + 
				GdPtCoef2(Rule, Grid) * Node(Node2) + 
				GdPtCoef3(Rule, Grid) * Node(Node3);
		}
}

//Set the integral coefficient**************************************************************************
void MoM_Geometry::SetCoef() 
{
	TinyVector <real_data ,16> AxisU,AxisV,Weight;
	for (int Rule=0; Rule<2; Rule++)
	{
		switch (Rule) 
		{
		case 0:
			AxisU = 0.0;
			AxisV = 0.0;
			Weight = 1.0;
			FillGrid(Rule, 1, AxisU, AxisV, Weight);
			break;
		case 1:
			AxisU = 0.0, 0.5, -0.25, -0.25;
			AxisV = 0.0, 0.0, 0.25 * SQRT3, -0.25 * SQRT3;
			Weight = 1.0/4.0, 1.0/4.0, 1.0/4.0, 1.0/4.0;
			FillGrid(Rule, 4, AxisU, AxisV, Weight);
			break;
		}
	}
}

//set the grid points' coefficients and weights*********************************************************
void MoM_Geometry::FillGrid(int Rule, int IntegralNum, const TinyVector<real_data ,16>& AxisU, 
	const TinyVector<real_data ,16>& AxisV, const TinyVector<real_data ,16>& Weight)
{
	Grid_Num(Rule) = IntegralNum;
	for (int IntegralPt=0; IntegralPt<IntegralNum; IntegralPt++)
	{
		GdPtCoef1(Rule, IntegralPt) = (1.0 - AxisU(IntegralPt) + SQRT3 * AxisV(IntegralPt)) / 3.0;  
		GdPtCoef2(Rule, IntegralPt) = (1.0 - AxisU(IntegralPt) - SQRT3 * AxisV(IntegralPt)) / 3.0;
		GdPtCoef3(Rule, IntegralPt) = (1.0 + AxisU(IntegralPt) * 2.0) / 3.0;
		GdPtWt(Rule, IntegralPt) = Weight(IntegralPt);
	}
}

